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Abstract 

Using molecular dynamic simulations we study a family of continuous core-softened potentials 
consisting of a hard core, a shoulder at closest distances and an attractive well at further distance. 
The repulsive shoulder and the well distances represent two length scales. We show that if the 
first scale, the shoulder, is repulsive or has a small well, the potential has a region in the pressure- 
temperature phase diagram with density, diffusion and structural anomalies. However, if the closest 
scale becomes a deep attractive well the regions in the pressure-temperature phase diagram where 
the three anomalies are present shrink and disappear. This result enables us to predict by the 
shape of the core-softened potential if anomalies would or would not be present. 

PACS numbers: 64.70.Pf, 82.70.Dd, 83.10.Rs, 61.20.Ja 



1 



I. INTRODUCTION 



Most liquids contract upon cooling. This is not the case of water, a liquid where the 
specific volume at ambient pressure starts to increase when cooled below T ^ 4°C— . Besides, 
in a certain range of pressures, also exhibits an anomalous increase of compressibility and 
specific heat upon cooling^i^i^. Experiments for Te,^ Ga, Bi,- S,-!^ and GeisTess,— and 
simulations for silica,— li^ii^i^ silicort^ and BeF2,— show, as well, the same density anomaly. 

Water also has dynamic anomalies. Experiments show that the diffusion constant, D, 
increases on compression at low temperature, T, up to a maximum D^.^^{T) at p = PDmax(^)- 
The behavior of normal liquids, with D decreasing on compression, is restored in water only 
at high p, e.g. for p > pumax ~ 1-1 kbar at lO^'C-'^ Numerical simulations for SPC/E water- 
recover the experimental results and show that the anomalous behavior of D extends to the 
metastable liquid phase of water at negative pressure, a region that is difficult to access 
for experiments.— li^ii^i^S In this region the diffusivity D decreases for decreasing p until it 
reaches a minimum value Z)min(T) at some pressure pDmin{T), and the normal behavior, 
with D increasing for decreasing p, is reestablished only for p < PDmin{T)^^^^^^- Besides 
water, silicap^i2^ and silicon^^ also exhibit a diffusion anomalous region. 

It was proposed a few years ago that these anomalies are related to a second critical point 
between two liquid phases, a low density liquid (LDL) and a high density liquid (HDL)^. 
This critical point was discovered by computer simulations. This work suggests that this 
critical point is located at the supercooled region beyond the fine of homogeneous nucleation 
and thus cannot be experimentally measured. Even with this limitation, this hypothesis has 
been supported by indirect experimental results^^i^^. 

In order to describe the anomalies present in water and in other liquids, isotropic models 
has been used as the simplest framework to understand the physics of the liquid-liquid 
phase transition and liquid state anomalies. From the desire of constructing a simple two- 
body isotropic potential capable of describing the complicated behavior present in water-like 
molecules, a number of models in which single component systems of particles interact via 
core-softened potentials'^ have been proposed. They possess a repulsive core that exhibits a 
region of softening where the slope changes dramatically. This region can be a shoulder or a 
ramp^^^^^>23^^^^^^J04i^4344^^^48^_ These models exhibit density, diffusion 

and structural anomalies, but depending on the specific shape of the potential, the anomahes 
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might be hidden in the metastable and unstable phases^. The relation between the specific 
shape of the core-softened potential and the presence or not of the anomalies is still missing. 

How the specific shape of a core-softened potential affects the location of the anomalies 
and the critical points? In order to answer to this question in this paper we analyze a family 
of continuous core-softened potentials that exhibit two length scales, a shoulder followed 
by an attractive well. When the shoulder is purely repulsive, this core-softened potential 
represents the effective pair interaction between two neighbors tetramers^*^ and the den- 
sity, the diffusion and the structural anomalies are present^^. If the shoulder has a deep 
attractive well, this potential it is related to the effective interaction potential between two 
water molecules obtained from the ST4^^ or TIP5P— models for water. In this case the 
effective potential is derived from the oxygen-oxygen radial distributions function, solving 
the Ornstein-Zernike equation by using an integral equation method^i^. The resulting po- 
tential has a shoulder with a deep attractive well at closest distance and a second attractive 
well with lower energy at furthest distance. The detailed depth of the softening region de- 
pends on the approximations employed. This potential leads, as we are going see in this 
paper, to systems in which the anomalies are in the unstable region of the phase diagram 
while in the full ST4 and TIP5P systems the anomalies can be observed. It is important, 
therefore, to understand what is lost when one goes from the specific anisotropic ST4 and 
TIP5P potentials to the isotropic spherical symmetric case. 

So, in this paper we study what happens with the region in the pressure-temperature 
phase diagram where the anomalies are located as the potential changes from a repulsive 
shoulder to a very deep well. Our results will shade some light not only in the use of spherical 
symmetric approximations of asymmetric potentials but also will help to design potentials 
for new systems with anomalies. 

The paper is organized as follows. In sec. II the family of potentials is introduced and 
its link with the derivation the framework of the integral equations is presented. In sec. Ill 
these potentials are tested for presence density, diffusion and structural anomalies, and for 
the the existence of two liquid phases and a critical points by molecular dynamic simulations. 
Conclusions are presented in sec. IV. 
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II. THE MODEL 



We study a system of particles, with diameter a, where the pair interaction is described 
by a family of continuous potentials given by 




* 



r 

FIG. 1: Interaction potential obtained by changing parameters hi in Eq. ([T]). 

The first term is a Lennard- Jones potential-like and the second one is composed by four 
Gaussian, each one centered in cj. This potential can represent a whole family of inter- 
molecular interactions, depending of the choice of the parameters a,b,a, {hj,Cj,Wj}, with 
j = 1, . . . , 4. The parameters are chosen in order to obtain a two length scale potential^. 

Modifying hi in the Eq. [Hallow us to change the depth of the hard-cor e well, as illustrated 
in Fig. [H Here we use four different values for hi and they are expressed as a multiple of a 
reference value h^^^^ as shown in the Table [H For all the four cases the values of a, b, {cj, wj} 
with j = 1, . . . , 4 and h'^'^^ are given in the Table [Tll The depth of the region of softening of 
the potentials illustrated in the Fig. [T] where chosen so that the potential B is the shallow 
shoulder-like potential similar to the one studied by de Oliveira et al.— that exhibits the 
anomalies, while for the potential D the region of softening has the same depth as the 
potential obtained by using the oxygen-oxygen radial distribution function for the ST4 
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TABLE I: Parameters hi for potentials A, B, C and D. 



Potential 


Value oi hi 


A 


0.25 h\'f 


B 


0.50 /i^^ 


C 


Q.15h\'f 


D 


1.00 /i^^ 



model^. For comparison we also analyzed two other cases: potential A with a ramp-like 
shoulder and potential C, with a very shallow shoulder. 

The properties of the system were obtained by NVT molecular dynamics using Nose- 
Hoover heat-bath with coupling parameter Q = 2. The system is characterized by 500 
particles in a cubic box with periodic boundary conditions, interacting with the intermolec- 
ular potential described above. All physical quantities are expressed in reduced units and 
defined as 

t(m/7e)^/2 



P* 

P* 
D* 



a 



7e 
pa 

7e 

pa' 

Dm 

e'ja'^ 



where 7 = 50. Standard periodic boundary conditions together with predictor-corrector 
algorithm were used to integrate the equations of motion with a time step At* = 0.002 and 
potential cut off radius r* = 3.5. The initial configuration is set on solid or liquid state 
and, in both cases, the equilibrium state was reached after t*^ = 1000. From this time on 
the physical quantities were stored in intervals of Atjj. = 1 during t*^^ = 1000. The system 
is uncorrelated after = 10, from the velocity auto-correlation function. 50 descorrelated 
samples were used to get the average of the physical quantities. The thermodynamic stability 
of the system was checked analyzing the dependence of pressure on density, by the behavior 
of the energy and also by visual analysis of the final structure, searching for cavitation. 



5 



TABLE II: Parameters for potentials A, B, C and D. 



Parameter 


Value 


Parameter 


Value 


a 


9.056 


Wl 


0.253 


b 


4.044 




1.767 


c 


0.006 


W3 


2.363 


d 


4.218 




0.614 


Cl 


2.849 




-1.137 


C2 


1.514 


h2 


3.626 


C3 


4.569 


hs 


-0.451 


C4 


5.518 


/14 


0.230 



III. RESULTS 

Pressure- Temperature Phase Diagram 

First, we are going to show the effects of the shoulder depth in the presence or not of 
the thermodynamic anomalies and the location in the pressure-temperature phase diagram 
of the different phases. Fig. [2] illustrates the pressure-temperature phase diagram of the four 
studied cases. The system has at high temperatures a fluid phase and a gas phase (not 
shown). These two phases coexist at a first order line that ends at a critical point (see 
Table IIIII for the pressure and the temperature values). At low temperatures and high 
pressures there are two liquid phases coexisting at a first order line ending at a second 
critical point (see Table IIVI for the pressure and the temperature values) that is identified 
in the graph by the region where isochores cross. 

In the Fig. [2] at low temperatures and low pressures the dotted line separates the fluid 
phase from the amorphous region where the diffusion becomes zero. For the potential A, 
the amorphous region is located in a pressure range —0.61 ^ p* ^ 3.40, for B case this 
region is located in the range —0.75 ^ p* ^ 0.40 and for C case it is located in the range 
—0.75 > > 0.40. The potential D does not has a stable amorphous phase. Hence, as the 
shoulder becomes deeper the amorphous phase shrinks and moves to a lower pressure range. 

6 



At low temperatures and high pressures two hquid phases are present. As the shoulder 
becomes deeper the liquid-liquid coexistence line slides down to lower pressures and it goes 
to higher temperatures. This indicates that the deeper the shoulder the liquid-liquid phase 
transition stays stable for higher temperatures. Therefore, even thought this transition only 
exists if the attractive part of the potential is present (the second length scale), the stability 
of the liquid phases is determined by the depth of the shoulder (the first length scale) . 



P 




FIG. 2: Pressure-temperature phase diagram for cases A, B, C and D. The thin solid lines are the 
isochores 0.30 < p* < 0.65. The liquid-liquid critical point is shown as a dot, the temperature of 

maximum density is a solid thick line, the diffusion extrema is the dashed hne and the structural 
extrema is the dashed-dotted line. The dotted line indicates the limit between the fluid and the 
amorphous regions. 
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TABLE III: Critical point location for potentials A, B, C and D. 



Potential 




Pel 


A 


1.93 


0.072 


B 


1.98 


0.078 


C 


2.02 


0.080 


D 


2.15 


0.094 



TABLE IV: Second critical point location for potentials A, B, C and D. 



Potential 




Pc2 


A 


0.35 


3.44 


B 


0.48 


1.86 


C 


0.57 


0.49 


D 


0.81 


-0.33 



Thermodynamics anomaly 

The Fig. [2] also shows the isochores 0.30 < p* < 0.65 represented by thin solid lines. The 
temperature of maximum density at constant pressure coincides with the minimum pressure 
on isochores, (^) = 0. From the equation 



dV 
df 



dp\ /'dV\ 



(2) 



is possible to see that, for a fixed density, a minimum in the pressure as a function of 
temperature represents a maximum in the density as a function of temperature, named 
temperature of maximum density (TMD) given by (1^) = 0. The TMD is the boundary 
of the region of thermodynamic anomaly, where a decrease in the temperature at constant 
pressure implies an anomalous increase in the density and therefore an anomalous behavior 
of density (similar to what happens in water). Fig. [2] shows the TMD as a thick solid line. 
For the potentials A, B and C the TMD is present but for potential D no TMD is observed. 

Similarly to what happens with the location of amorphous region and of the second crit- 
ical, as the shoulder becomes deeper, the region in the pressure-temperature phase diagram 
delimited by the TMD goes to lower pressures, shrinks and disappears for the case the 
potential with the deepest shoulder. As the region delimited by the TMD shrinks, it also 
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FIG. 3: Location of the critical points on pressure-temperature phase diagram for cases A, B, C 
and D. 

goes to lower temperatures. For the potential C the temperature range of the TMD is lower 
than the liquid-liquid critical point. The thermodynamic parameters that limits the TMD 
in phase diagram are shown in Table |Vl where pi represents the pressure of the point with 
the lowest pressure, ph the pressure of the point with highest pressure and Pm the pressure 
of the point with the highest temperature. 

The link between the depth of the shoulder and the presence or not of the TMD goes as 
follows. The TMD is related to the presence of large regions in the system in which particles 
are in two preferential distances represented by the first scale and the second scale in our 
potential iEi^iMi^. While for normal liquids as the temperature is increased the percentage 
of particles at closest scales decreases (see case D in the Fig.H]), for the anomalous liquid (see 
cases A, B and C in the Fig. Hj) there are a region in the pressure-temperature phase diagram 
where as the temperature is increased the percentage of particles at the closest distance 
increases. Particles move from the second to the first scale. In the first case, the decrease 
of particles in the first scale leads to a decrease of density with increase of temperature, 
behavior expected for normal liquids. In the second case, the increase of particles in the 
first scale leads to an increase of density with temperature what characterizes the anomalous 
region. Notice that as the temperature is increased particles move from the second scale 
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cases 


Pi 


Pm 


Ph 


P* 


0.47 


0.52 


0.57 


A T* 


0.71 


0.85 


0.73 


P* 


1.50 


2.50 


3.30 


P* 


0.46 


0.50 


0.54 


B T* 


0.67 


0.76 


0.63 


p* 


0.90 


1.40 


1.80 


P* 


0.40 


0.42 


0.43 


C T* 


0.44 


0.54 


0.52 


* 

P 


0.15 


0.29 


0.36 



TABLE V: Limits values for density (p*), temperature (T*) and pressure (p*) of the thermody- 
namics anomahes on pressure-temperature diagram. Where pi is a lower limit (lesser pressure), Pm 
is a inflection point of the anomaly (higher temperature) and ph is a higher limit of the anomalies 
(higher pressure). 

to the first scale. The anomaly is, therefore, related with the possibility of having particles 
moving from one scale to the other as temperature is changed what becomes quite difficult 
if the depth of the shoulder well becomes too deep. 



Diffusion anomaly 

Now we are going to test the effect the shoulder depth has in the location of the diffusion 
anomaly in the pressure temperature phase diagram. The diffusion coefficient is obtained 
from the expression: 

D = lim + -^ifa)]^)*o (^) 

where rj{t) are the coordinates of particle j at time t, and (■ ■ ■ )tQ denotes an average over 
all particles and over all to- 

Fig. [5] shows the behavior of the dimensionless translational diffusion coefficient, D* ^ as 
function of the dimensionless density, p*, at constant temperature for the four cases. The 
solid lines are a polynomial fits to the data obtained by simulation (the dots in the Fig. [5]). 
For normal liquids, the diffusion at constant temperature increases with the decrease of 
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FIG. 4: Radius distribution as a function of the distance for the four potentials. In the cases 
B and C the first peak of g(r*) grows in increasing temperature, while the second peak decreases. 
For the potential D all the peaks decreases with temperature. 

the density. For the potentials A, B and C the diffusion has a region in the pressure- 
temperature phase diagram where the diffusion increases with density what represents a 
diffusion anomalous region. In the Fig. [5] one dashed line joints the points of the density (or 
pressure) of minimum diffusion for different temperatures and another dashed line links the 
points of density (or pressure) of maximum diffusion for different temperatures. 

Similarly to what happens with the location of the TMD, as the shoulder becomes deeper, 
the region in the pressure-temperature phase diagram delimited by the extrema of the dif- 
fusion goes to lower pressures, shrinks and disappears for the case the potential with the 



11 




0,4 0.5 0.6 0.7 

P 



0.4 0.5 0.6 0.7 

P 




0.4 0.5 0.6 0.7 

* 

P 



0.8 



FIG. 5: Diffusion coefficient as a function of density. Tlie dots are tlie simulational data and the 
solid lines are polynomial fits. The dashed lines connect the densities of minima and maxima 
diffusivity that limit the diffusion anomalous region. 

deepest shoulder. 

Fig. [2] shows the location at the pressure-temperature phase diagram of the pressure of 
maximum and minimum diffusion as double dot dashed lines (the dashed lines in the Fig. [5]). 
In the Fig. [2] we show that in the pressure-temperature phase diagram the region where the 
dynamic anomaly occurs englobes the region where the thermodynamic anomaly is present. 
This hierarchy between the anomalies is observed in a number of models^^'i^'^ and in the 
water.- 

The link between the depth of the shoulder and the presence or not of the region of 
diffusion extrema goes as follows. The presence of the diffusion anomaly is related to having 
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the quantity S2 > 0A2M^ where 

S2 = 



dS2 

dlnp 



(4) 



S2 - 2ttp 



(jQiv] 

In q(r) — - — r'^dr 
^ ' dp 



where 



S2 = -27TP j [g{r) \ng{r) - g{r) + 1] r'^dr , (5) 

is the excess entropy. Fig. [6] illustrates the behavior of the radial distribution function for 
fixed temperature as the density varies. For the case A the hig{r) is negative and dg{r) / dp is 
positive for the first scale, while for the second scale the hig{r) is positive and the dg{r)/dp 
is negative. As a result the second parcel in Eq. ([5]) is positive a requirement for having 
S2 > 0.42 since S2 is negative^. For case D, also shown in Fig. [6l the In g{r) is positive 
and huge and dg{r)/dp is positive what leads to a second parcel in Eq. ([5]) that is negative 
what do not fulfill the requirement S2 > 0.42. If the shoulder well is too deep the particle 
is unable to go from one scale to the other and the density anomalous behavior does not 
h 



p* = 0.400 

__ p* = 0.430 
• p* = 0.460 
p* = 0.570 




FIG. 6: Radial distribution for cases A and D as a function of r* for various densities. In the case 
A the temperature is fixed T* = 0.90 while in the case D the temperature is T* = 1.10. 



Structural anomaly 



Finally we are going to test the effect the shoulder depth has in the location in the 
pressure-temperature phase diagram of the structural anomalous region. 
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FIG. 7: The translational order parameter as a function of density for fixed temperatures: T* = 
1.10, 1.00, 0.90, 0.80, 0.70 and 0.60 (from top to bottom). The dot-dashed lines locate the density 
of maxima e minima t* . 

The translational order parameter is defined sisM^^ 

t= [^^m)~MdC (6) 
Jo 

where ^ = rp5 is the distance r in units of the mean interparticle separation p^^, is the 
cutoff distance set to half of the simulation box times^ P~K diO is the radial distribution 
function proportional to the probability of finding a particle at a distance ^ from a referent 
particle. The translational order parameter measure how structured is the system. For an 
ideal gas g = I and t = 0, and the case of crystal phase g ^ 1 over long distances and t is 
large. Therefore for normal fluids t increases with the increase of the density. 

Fig. [7] shows the translational order parameter as a function of density for fixed temper- 
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ature. The dots represent the simulation data and the sohd hne the polynomial fit to the 
data. For the potentials A, B and C there are a region of densities in which the translational 
parameter decreases as the density increases. A dotted-dashed line illustrates the region of 
local maximum of t* and minimum of t* limiting the anomalous region. For the potential 
D, t* increases with the density. No anomalous behavior is observed. 

Fig. [2] shows the structural anomaly for cases A, B and C, as dashed-pointed (purple) 
lines. It is observed that the region of structural anomaly embraces both dynamic and 
thermodynamic anomalies. Similarly to other anomalies the effect of increase the depth of 
the repulsive shoulder is to narrow the anomalies asymmetrical. The branch of anomaly 
in pressures near to liquid-liquid critical point is most feeling to the effect of the shoulder 
compared with the branch obtained in low pressures. However, the hierarchy of the anomalies 
is maintained, the change in the repulsive shoulder does not affect it. 

IV. CONCLUSIONS 

In this paper we studied a family of potentials characterized by two length scales: a 
repulsive shoulder and an attractive well. We analyzed the effect in the location in the 
pressure-temperature phase diagram of the density, diffusion and structural anomalies of 
making this repulsive shoulder a deep well. We found that the anomalies shrink and dis- 
appear as the well becomes deeper. This indicates that an important mechanism for the 
anomalies is the possibility of particles in the furthest length scale to move to the closest 
length scale. As the shoulder well becomes deeper particles becomes localized in the closest 
scale and the mobility between the two scales decreases. 

We find that in the cases of potentials A, B and C the thermodynamic, dynamic an 
structural anomalies are present and that the region of structural anomaly embraces the 
dynamic and thermodynamic anomaly in pressure-temperature phase diagram. This implies 
that the hierarchy of the anomalies is preserved independent of the depth of the repulsive 
shoulder, however when the shoulder becomes deeper, the upper pressure lines of anomaly 
converge to a similar value in the pressure-temperature phase diagram. 

What is the connection between the studies potentials and the real system? Effective 
potentials for water has been derived based in the oxygen-oxygen radial distribution func- 
tion for the ST4^^ and TIP5P— models for water. In both cases the effective potential was 
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obtained from the g(r*) using the Ornstein-Zernike equation and integral equation approxi- 
mations. The potential resulting are the case D in the Fig. [1] in the case of ST4 and for the 
TIP5P model a potential that exhibits a deep shoulder similar to the case D. Consequently 
the approximation washes out the anomalies present in both ST4 and T1P5P. In the case 
of the TIP5P it was shown that if instead of deep shoulder a smooth shoulder like the one 
present in the ramp potential would be used, the anomalies not only would be present but 
would be located in the same region of pressure and temperature of the TIP5P potential. 

In resume, similarly to other previous studie8^i^>^i^, a directional interaction potential 
is not a fundamental ingredient to have a water-like anomalies. Two scales isotropic potential 
also reproduce this anomalies if the shoulder closest scale would not be too deep. 
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